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ABSTRACT 

We present the results of an ensemble of SPH simulations that follow the evolution 
of prestellar cores for 0.2 Myr. All the cores have the same mass, and start with the 
same radius, density profile, thermal and turbulent energy. Our purpose is to explore 
the consequences of varying the fraction of turbulent energy, that is solenoidal, 
as opposed to compressive; specifically we consider (igoL = 2/3, 1/3, 1/9 and 0. For 

each value of ^sol, we follow ten different realisations of the turbulent velocity field, in 
order also to have a measure of the stochastic variance blurring any systematic trends. 
With low ^sol(< 1/3) filament fragmentation dominates and delivers relatively high 
mass stars. Conversely, with high values of (^sol(> 1/3) disc fragmentation dominates 
and delivers relatively low mass stars. There are no discernible systematic trends in 
the multiplicity statistics obtained with different (igoL- 

Key words: 


I INTRODUCTION 

Understanding the origin of the stellar initial mass function 
(IMF) (e.g 


Kroupa 2001 Chabrier 2003 


2005) is one of 
the main unresolved problems in star formation. Because 
the physics regulating star formation (i.e. self-gravity, hy¬ 
drodynamics, radiative transfer, magnetism, etc.) is highly 
non-linear, general cases can only be studied with numeri¬ 
cal simulations. An intrinsic feature of the initial conditions 
for such simulations is the imposed turbulent velocity field, 
and since this velocity field is the source of the density fluc¬ 
tuations that spawn protostars, it is important to be clear 
about how it is defined. 

Simulations of star formation usually follow one of two 
approaches. The first approach involves the simulation of 
isolated prestellar cores, i.e. the small {R ~ 0.1 pc), dense 
(p > 3 X 10“^° g cm“^) clumps of gas with subsonic or mildly 
transsonic internal turbulence, in which individual stars 
or small sub-clusters form (e.g 
Bate & Bonnelll 


Whitworth |2004[ |Delgado-Donate, Clarke fc Bate||2004| 


2000 


Matsumoto fc Hanawa|200^ Goodwin 




Horton, 


Delgado-Donate et al. 2 

i004| Goodwin, Whitworth & Ward- 

Thompson||2004| 

2006 

Walch et al.||2009[ Federrath et al. 

20 10a[| Walch et a 

l.|201C 

Girichidis et al.|2011||Walch, Whit- 


worth Girichidis|2012 ). This approach has the advantage 

that individual simulations can be performed at high res- 


E-mail: oliver.lomax@astro.cf.ac.uk 


olution with modest computational resource. Gonsequently 
good statistics can be obtained by performing multiple re¬ 
alisations. However, it is still important to ensure that the 
initial conditions mimic reality as closely as is possible (e.g. 
Lomax et al.|2014 ). 

The second approach involves the simulation of the 
much larger, less dense and more massive molecular clouds 
with highly supersonic internal turbulence, in which prestel¬ 
lar cores form (e.g. Klessen, Heitsch fc Mac Low|2QQQ) Bate 


2009 2012 Federrath Klessen|2012 Bate|2014 ). This ap¬ 

proach has the advantage that the evolution includes both 
the formation of cores, and interactions between them, but 
it is not always feasible to perform more than one realisa¬ 
tion, and it remains to be seen if the initial conditions on 
this scale are critical. 

In this paper we consider the mildly transsonic turbu¬ 
lent velocity fields used to initiate simulations of individual 
prestellar cores. Specifically, we explore the effect of chang¬ 
ing the fraction of turbulent energy, (5 sol, that is solenoidal 
as opposed to being compressive. Techniques for measuring 
the ratio of solenoidal to total turbulent energy have re¬ 
cently been developed ( Brunt fc Federra'th||2Q14 ), but have 
not yet been widely applied. The values of Jsol invoked in 
numerical simulations are seldom explicitly justified. The 
most common choices are Ssol = 2/3 (thermal mixture, e.g. 
Walch, Whitworth fc Girichidis|[2012| [Lomax et "ar 2014) 


and 6sol = 1 (purely solenoidal, e.g. |Bate|2009||2012|(2014| | 

Previous numerical work on highly supersonic turbu- 
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lence on molecular cloud scales (e.g. Federrath & Klessen 


2012 ) indicates that compressive turbulence can deliver star 


formation rates up to ten times higher than solenoidal tur¬ 
bulence. Numerical simulations of very massive cores with 
M = 100M© and R = 0.1 pc (Girichidis et al. 2011) also 
show that compressive turbulence accelerates the onset of 
star formation relative to solenoidal turbulence. Here we 
study this issue on the much smaller scale of the prestel- 
lar cores typically seen in nearby star forming regions (e.g. 
Ophiuchus), where the turbulence is a lot less vigorous. 

In ^we detail the initial conditions used for the sim¬ 
ulations. In E we describe the numerical method and the 
constitutive physics. In ^we present the results, and in ^ 
we summarise our conclusions. 


2 INITIAL CONDITIONS 

All the simulations presented here start with a spherical 
core having total mass M = 3 M©, radius R = 3000 an and 
non-thermal velocity dispersion ctnt = 0.44kms“^.These 
values are similar to those of the SMI core in the Oph- 
A clump within the LI688 (Ophiuchus) cloucj^ The initial 
background temperature is set to 10 K. 


2.1 Modified random turbulent velocity field 


2.1.1 Standard random turbulent velocity field 

Following Lomax et al. ( 2014| , each core initially has a tur¬ 
bulent velocity_field with power spectrum Pk oc k~^ (Burg¬ 


ers turbulenc Cl where k — SttR/X is the wavenumber of a 
velocity mode having wavelength A. In three dimensions, 
each velocity mode is characterised by (i) a wavevector 
k = (/ci,A: 2 ,fe); (ii) an amplitude 


a(fc) = , (2.1) 


where the Qn are random variates from a Gaussian distri¬ 
bution (mean /x = 0, standard deviation cr = 1); and (iii) a 
phase 


(UA 

¥>(fc) = 27t , (2.2) 

\UsJ 


where the Un are random variates from a uniform distribu¬ 
tion on the interval [0,1]. Non-zero amplitudes are given to 
wavevectors with integer components satisfying 


0 ^ /cs < kuAx 5 


^MAX ^ /C2 ^MAX 

if ks > 

0 < /C2 ^MAX 

if ks = 

^MAX kl /CmAX 

if ks > 

0 ki /Cmax 

if ks = 


0, 

0, (2.3) 

0 or /c 2 > 0, 

0 and /c 2 = 0. 


^ Using 1.3 mm dust continuum observations, [Motte, Andre &:| 
|Neri| (|1998| estimate that SMI has a mass of 3.2 M© and an az- 
imnthally av eraged full width at half-maximum of 3600 au. |Andre| 
|et al.| ( [2007| measure the velocity width of the N 2 H+ (1-0) line 
and estimate that the three-dimensional non-thermal velocity dis¬ 
persion is 0.45 km s“^. 


These wavevectors cover all frequencies up to the Nyquist 
frequency of a grid with (2 /cmax)^ uniformly spaced elements. 


2.1.2 Modifications to the longest wavelength modes 


The longest-wavelength velocity modes, those correspond¬ 
ing to the scale of the core, are then modified so that they 
generate radial excursions (either contraction or expansion) 
relative to the centre of the core, and rotation about the cen¬ 
tre of the core. This is achieved by revising the amplitudes 
of the modes ki = (1, 0,0), k^ — (0,1,0) and k^ — (0,0,1) 
to 


(g.\ 

a(fci) = Qe , a{k- 2 ,) 

\-gS 



and their phases to 



(2.4) 


(n/2\ 

^(ki) = (p{k2) = (fiks) = 71/2 . (2.5) 

v/v 


With this procedure, ^i, Q 2 and determine the amount 
of global compression (expansion) towards (away from) the 
centre of the core; and ^ 4 , Qb and Qq determine the amount 
of global rotation about the centre of the core. This adjust¬ 
ment is made because, if a newly-formed core is undergoing 
global contraction or expansion, these motions are likely to 
be focussed on the centre of the core. Similarly, if a newly- 
formed core is undergoing global rotation, these motions are 
likely to be around the centre of the core. All other veloc¬ 
ity modes, up to /cmax = 64 retain their random phases, as 
generated by Eqns. (2.1) and (2.2), and therefore represent 
internal random turbulence. 


2.1.3 Helmholtz decomposition 

Helmholtz’s theorem states that a vector field can be ex¬ 
pressed as the sum of a compressive (curl-free) vector field 
and a solenoidal (divergence-free) vector field. For a velocity 
mode with wavevector k and amplitude a{k), the longitu¬ 
dinal component of the amplitude, 

cii^{k) = k{k • a{k)) , (2.6) 

contributes to the compressive field Vc{x), and the trans¬ 
verse component of the amplitude, 

aT{k) = a{k) — ai^{k) , (2.7) 

contributes to the solenoidal field Vs{x). These components 
can be summed as wave amplitudes or in real space, 

a(fc) = aPk) + ax(fc), 
v{x) = Vcix) + Vsix), 


^ Strictly speaking, a power law exponent of -4 (Burgers turbu¬ 
lence) is only appropriate for highly supersonic turbulence (see 
|Federrath]|2013| ). At sonic and transonic speeds, the exponent 
is likely to be between -4 and -11/3 (Kolmogorov turbulence). 
However, with both exponents, the turbulent energy is strongly 
concentrated at the longest wavelengths, and so the precise choice 
of exponent is not critical. 
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to retrieve the total values. Note that ai^(k) is the compo¬ 
nent of a{k) parallel to k and aT{k) is the perpendicular 
component. In three dimensions - and assuming that the 
d{k) are distributed isotropically - Vs(x) has on average 
twice the kinetic energy of Vc{x). This is because trans¬ 
verse waves have two degrees of freedom whereas longitudi¬ 


nal waves only have one (see Federrath, Klessen & Schmidt 


2008). Helmholtz decomposion has also been used in other 
astrophysical simulations (e.g. Schmidt et al. 2009| Feder- 


rath et al.|2010b Girichidis et al.|2011). 


For an arbitrary velocity field v{x), we can alter the 
fraction of solenoidal kinetic energy by decomposing and 
then reconstituting the amplitudes of each velocity mode. 
Here, we define five sets of modified amplitudes: 


ai(fc) = ttT^k) , 

a^ik) — ai^{k) + ar^ik) , 

as{k) = 2ai^{k) + aT{k) , (2.9) 

a^ik) = 4aL(fc) + aT(k) , 

a^ik) = ai^{k) . 


The average fractions of kinetic energy in solenoidal 
modes for the corresponding velocity fields, vi(x), 
V 2 {x), vsix), V 4 :{x) and v^ix), are respectively Ssol = 
1, 2/3, 1/3, 1/9 and 0. 


2 . 1.4 Particle velocities 

Having defined all the velocity modes, we use the fast 
Fourier transform library FFTW (Frigo & Johnson 2005) 
with /cmax = 64 to compute a gridded velocity field v{x) on 
—2R < xi,X 2 ,X 3 < +2i?, and the velocities from the cen¬ 
tral eighth of the volume are then mapped onto the SPH 
particles. 


2.2 Density profile 


Many observations (e.g. |Alves, Lada & Lada| 

2001 

Har- 

vey et al.|2001[ |Kirk, Ward-Thompson & Andre 

2005 

|Lada 

et al.|2008| Roy et al.|2014|) suggest that the critical Bonnor- 


Ebert sphere provides a good fit to the column-density pro¬ 
file of a prestellar core, even if the core is not in hydro¬ 
static equilibrium. We therefore set the core density profile 
to p(^) = pce~'^^^\ where pc is the central density, '0(^) is 
the Isothermal Function and ^ is the dimensionless radius, 
i.e. = 6.451(r/3000au). 


2.3 Parameter space 

Using the procedures described in Sections |2TT] and [2T^ 
we generate ten different initial velocity fields, by invoking 
ten different random seeds, 

Tseed = 1, 2 , 3, 4, 5, 6, 7, 8, 9 and 10. (2.10) 

Then, using the procedures described in Section |2.1.3[ we 
convert each of these velocity fields into five velocity fields 
with different fractions of solenoidal kinetic energy 

JsoL = 1, 2/3, 1/3, 1/9 and 0. (2.11) 

We therefore have a total of 50 initial velocity fields, cor¬ 
responding to all possible combinations of the ten different 
random seeds, Xseed, and the five different fractions, Jsol- 



Figure 1. False-colour column-density images on the central 
820 au by 820 au of the (x, y)-plane, from the simulation with 
Xseed = 3 and Jsol = 2/3, at times t = 0.85, 0.90, 0.95 and 1.00 x 
10"^ yrs. The colour scale gives the logarithmic column density in 
units of gcm“^. Sink particles are represented by black dots. This 
is an example of filament fragmentation, where the filaments serve 
to deliver matter from the periphery of the core into the centre. 
Further evolution of this case is shown in Fig. |^. 


3 NUMERICAL METHOD 


Core evolution is simulated using the seren V/i-SPH code 
(Hubber et al. 1 |MT] i, with ?7 = 1.2 (so a particle typically 
has 57 neighbours). Gravitational forces are computed us¬ 
ing a tree, and the Morris & Monaghan (1997) formulation 


of time dependent artificial viscosity is invoked. In all sim¬ 
ulations, the SPH particles have mass ttisph = 10“^ M©, so 
the opacity limit (~ 3 x 10“^ M©) is resolved with ~300 
particles. Gravitationally bound regions with density higher 
than psiNK = 10“^gcm“^ are replaced with sink particles 
( Hubber, Walch Whitworthj2013 ). Sink particles have ra¬ 
dius rsiNK — 0.2 au, corresponding to the smoothing length 
of an SPH particle with density equal to Psink- The equa¬ 
tion of state and the energy equation are treated with the 


algorithm described in Stamatellos et al. (2007). 


Radiative feedback from sinks is also included. Each 
sink has a variable luminosity which follows the episodic ac- 


cretion model described in|Stamatellos, Whitworth & Hub- 

ber ( 

201 1|) (also used in Stamatellos, Whitworth & Hubber 

2012 

Lomax et al. 

2014). In this model, highly luminous. 


short-lived accretion episodes are separated by ~ 10 yrs 
of low-luminosity quiescent accretion (during which matter 
collects in the inner accretion disc until it is hot enough 
to become thermally ionised and couple to the magnetic 
field; then the Magneto-Rotational Instability delivers effi¬ 
cient outward angular momentum transport and the matter 
is dumped onto the star). 
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Figure 2. False-colour column-density images on the central 
820 au by 820 au of the (x, y)-plane, from the simulation with 
^SEED = 3 and (5 sol = 0, at times t = 0.85,0.90,0.95 and 1.00 x 
10"^ yrs. The colour scale gives the logarithmic column density in 
units of gcm“^. Sink particles are represented by black dots. This 
is an example of filament fragmentation, where the individual fil¬ 
aments fragment independently to produce an ensemble of stars. 
Further evolution of this case is shown in Fig. |5(ay| 



Figure 3. False-colour column-density images on the central 
820 au by 820 au of the (a:, ^)-plane, from the simulation with 
Tseed — 3 and (5 sol — 2/3, at times t = 1.00,1.05,1.10 and 1.15 x 
10^ yrs. The colour scale gives the logarithmic column density in 
units of gcm“^. Sink particles are represented by black dots. This 
is an example of disc fragmentation. Further evolution of this case 
is shown in Fig. |5(d)] . 


RESULTS 


Each core is evolved for 0.2 Myr. This is roughly the pre¬ 


dicted time between core-core collisions in Ophiuchus (An- 
dre et al.|2QQ7 ). The simulations do not include any mechan¬ 


ical feedback from sinks (e.g. outflows), so the star formation 
efficiency is high, roughly 0.8 < rj < 1.0. Previous work (e.g. 
Matzner McKee|2000| Federrath et al.|2014 ) demonstrates 


that outflows and jets can reduce star formation efficiency 
signficantly. 


4.1 Modes of fragmentation 

As a core collapses, turbulence and self-gravity organise the 
matter into filaments. These filaments usually behave in one 
of two ways: (i) they feed material from the periphery of 
the core into its centre, forming a central star, or (ii) they 
fragment independently to form multiple stars, which then 
congregate in a small cluster near the centre of the core. An 
example of case (i) is shown in Fig. [^and an example of case 
(ii) is shown in Fig. In the sequel we refer to this initial 
mode of star formation as filament fragmentation. 

Once filament fragmentation has occurred, the remain¬ 
ing matter in the core envelope tries to accrete onto the star 
or cluster near the centre. If this matter has sufficient angu¬ 
lar momentum, it forms circumstellar or circumsystem discs. 


Discs which are sufficiently 

massive ( 

Toomre||1964) and are 

able to cool sufficiently fast (Gammie 

2001 ) fragment to pro- 

duce additional stars (e.g. 

Stamatellos & Whitworth||2008 

2009a|b Stamatellos et al. 

2011). An example of this pro- 


cess is shown in Fig.[^ In the sequel we refer to this second 
mode of star formation as disc fragmentation. 


4.2 Influence of (5 sol on the dominant mode of 
fragment at ion 

Table lists the number of stars that form by filament frag¬ 
mentation, A^ee, and the number that form by disc fragmen¬ 
tation, iVoF The distinction between the two modes is made 
by inspecting the simulation frames by eye. This table high¬ 
lights the need to invoke multiple realisations with different 
random seeds, since, with a given ^sol, the results can vary 
dramatically with Xseed- For example, with (5 sol = 1/9, only 
one star forms when Xseed = 5, whereas twelve stars form 
when Xseed = 10. Tffis is because most of the turbulent en¬ 
ergy is in large-scale modes which are defined by only a few 
wavevectors, and so the outcome is very sensitive to the ran¬ 
dom amplitudes of these modes. 

Fig. E shows the early stages of star formation with 
Xseed = 3 and different values of 6sol- In Fig. |5(a)| where 
Ssoh = 0, seven stars form by filament fragmentation, and 
then a single star by disc fragmentation. In Fig. |5(e)'| where 
(5 sol = 1 , a single star forms by filament fragmentation, then 
thirteen by disc fragmentation. Figs. |5(b)[ |5(c) and 5(d)| 
show how the number of stars formed by filament fragmen¬ 
tation tends to decrease, and the number of stars formed 
by disc fragmentation to increase, with increasing Ssoh- This 
trend is seen more clearly in Fig. where the results are 
averaged over all values of Xseed- The average fraction of 
stars formed from a core by filament fragmentation decreases 
monotonically from ~ 0.9 when (5 sol = 0 to ~ 0.2 when 
dsoh = 1- Conversely, the average fraction of stars formed 
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Figure 4. The fraction of stars formed by filament fragmentation 
(red xs) and disc fragmentation (green Ds) for different values of 
(^soL, averaged over all values of Xseed- The error bars give the 
Poisson counting uncertainties. 


from a core by disc fragmentation increases monotonically 
from ~ 0.1 when (5 sol = 0 to ~ 0.8 when (5 sol = 1 • This 
is because predominantly compressive fields (low (5 sol) are 
characterised by shocks, and these are conducive to filament 
formation. Conversely, predominantly solenoidal fields (high 
(5 sol) are characterised by shearing motions, and these gen¬ 
erate the angular momentum required for the formation of 
discs. 

Table [T] also indicates that the total number of stars 
formed per core increases slightly with increasing (5 sol. When 
(^soL = 0, a core spawns on average ~5 stars; when (5 sol = 1 , 
a core spawns on average ~8 stars. 


4.3 Influence of (5 sol on the stellar mass 
distribution 

Fig. E shows the distribution of stellar masses formed with 
different values of (5 sol, integrated over ah values of Xseed- 
Fig. shows how the corresponding medians and interquar¬ 
tile ranges vary with Ssoh- We see that, as (5 sol is increased, 
the median decreases monotonically, from ~ 0.6 M© when 
(^SOL — 0 , to ~0.3Mo when (5sol = 1. For reference, Fig. 
also shows the Chabrier (2005; hereafter C05) IMF (dashed 
red curve), and Fig.shows the corresponding median and 
interquartile range (full and dashed horizontal red lines). 
However, we stress that we should not expect to reproduce 
the overall distribution of stellar masses observed in na¬ 
ture with simulated cores of a single mass, radius and non- 
thermal velocity dispersion, as treated here. We are simply 
seeking to establish what trends, if any, might derive from 
changing the mix of solenoidal and compressive modes in 
the imposed turbulent velocity held. 

The decrease in median mass that results from increas¬ 
ing SsoL can be attributed directly to the shift from hlament 
fragmentation when (5 sol is low to disc fragmentation when 
Ssoh is high. When Ssoh is low, compressive turbulent modes 
create hlaments, and these can be very effective at feeding 
matter from the periphery of the core, into the centre, where 
it forms a relatively massive star. Conversely, when Ssoh is 
high, solenoidal turbulent modes create discs around exist¬ 
ing stars, and these discs tend to fragment to produce large 
numbers of low-mass companion stars. 



log(M/M0) 


Figure 6. The black histograms show un-normalised stellar mass 
distributions, integrated over all values of Xseed, for different val¬ 
ues of (5 sol- The error bars give the Poisson counting uncertainties. 
The red dashed lines show the |Chabrier| ( [200^ IMF, scaled to the 
area of the (5 sol = 1 histogram. 



^SOL 


Figure 7. The black xs give the median stellar mass, and the 
vertical black bars give the corresponding interquartile range, for 
different values of (5 sol, integrated over all values of Xseed- The 
solid and dashed horizontal red lines give the median and IQR 
for the |Chabrier| ( [^05| IMF. 
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<5sol 

= 0 

<5sol 
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<5sol 
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2 

4 

9 

4 

2 

4 

4 

1 

8 

2 

9 

1 

8 

10 

4 

0 

4 

8 

2 

5 

2 

9 

1 

5 

Total 

49 ±7 

5±2 

40 ±6 

22 ±5 

27 ±5 

39 ±6 

28 ±5 

60 ±8 

14 ±4 

67 ±8 

Fraction 

0.91 ±0.04 

0.09 ±0.04 

0.65 ±0.06 

0.35 ±0.06 

0.41 ±0.06 

0.59 ±0.06 

0.32 ±0.05 

0.68 ±0.05 

0.17 ±0.04 

0.83 ±0.04 


Table 1. The number of sinks formed by filament fragmentation, A^ff, and by disc fragmentation, A^df, in each simulation. Column 
1 gives the random seed. Columns 2 & 3 give the number of sinks formed by filament fragmentation and by disc fragmentation when 
(^soL = 0. Similarly, columns 4&5, 6&7, 8&9 and 10 & 11 give the same quantities when, respectively, (5 sol = 1/9, 1/3, 2/3 and 
l.The second from last row gives total number of stars over all random seeds. The last row gives the relative fraction of stars formed via 
filament and disc fragmentation. 




(d) (5 sol = 2/3 (e) (5sol = 1 


Figure 5. False-colour column-density images on the central 820 an by 820 an of the (a:, y)-plane, from the simulations with Xseed = 3 
and different values of (5 sol, at times t = 1.00,1.25,1.50 and 1.75 x 10^ yrs. The colour scale gives the logarithmic column density in units 
of gcm“^. Sink particles are represented by black dots. 
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(^SOL — 0 (^SOL — ^ 


NqyS 

32 

34 

mf 

0.28 

0.32 

Amf 

0.08 

0.08 

pf 

0.63 

0.82 

Apf 

> 0.18 

> 0.20 


(^SOL — ^ (^SOL — ^ (^SOL — 1 


51 

61 

55 

0.20 

0.23 

0.29 

0.06 

0.05 

0.06 

0.29 

0.42 

0.49 

> 0.08 

> O.IO 

> 0.10 


Table 2. Multiplicity frequencies and pairing factors, for differ¬ 
ent values of (5 sol, integrated over all values of Xseed- Uncertainties 
on mf are given by Am/ = ^/mf (1 — mf)/{NsYs) where Nsys 
is the total number of systems. Uncertainties on pf are assumed 
to satisfy Apf > pf Amf /mf (see footnote Section [4.4.1| for dis¬ 
cussion). 


4.4 Multiplicity statistics 

4-4-^ Multiplicity frequency 

Table lists the multiplicity frequencies and pairing factors 
extracted from the simulations. The procedures used to de¬ 


termine these statistics are detailed by Lomax et al. (2015). 


The multiplicity frequency, m/, is the fraction of systems 
which is multiple. The pairing factor, pf, is the mean num¬ 
ber of orbits per system (see Reipurth fc ZinneckerlpT^ ). 
Thus, if S is the number of single systems, B the number of 
binaries, T, Q, etc. the numbers of triples, quadruples, etc., 
then 


mf 

pf 


B + T + Q + ... 

5' + R + T + Q + ... 
R + 2T + 3Q + ... 
5' + R + T + Q + ... 


(4.1) 

(4.2) 


When all the multiple systems in a population are binary, 
pf = mf. When higher multiples are present (i.e. triples, 
quadruples, etc), pf > mf. We see from the table that there 
is no discernible trend in mf with changing Ssol- Values of 
mf range between 0.2 and 0.3, which is roughly the same as 
for M-dwarf stars in the field (see Duchene fc Kraus||2013| 
and references therein). In all cases, pf > mf, due to the 
presence of hierarchical multiple systems (up to sextuples). 
The variation of pf is seemingly much more stochastic than 
that of mf. However, the stated uncertainties of pf in Table 
l^are only lower limits 


4.4-2 Orbital properties 

The number of multiple systems formed in the simulations 
is too small to allow us to consider in detail how the distri¬ 
butions of orbital properties vary with (5 sol. There is some 
tennuous evidence that, as (5 sol increases, the orbital eccen¬ 
tricities increase and the mass ratios decrease - in other 


^ mf depends only on two numbers, S and B T Q = 
(Nqys — S). In contrast, pf depends on S, B, T, Q, etc.; if we 
count up to sextuples, then six numbers. The higher the order of 
the system being counted, O, the smaller the number of systems 
and hence the higher the fractional Poisson uncertainty. However, 
these systems are given higher weight in the numerator of pf, i.e. 

= (O — I), so their individual uncertainties are compounded 
in a way that they are not in mf. 


CTa 1 e 

0.7 ±0.1 I.OTO.I 1.6 ±0.2 1.6 ±0.2 


Table 3. Fitted multiplicity parameters integrated over all sim¬ 
ulations (i.e. all (5 sol and Xseed)- yta and da are the mean and 
standard deviation of log^^ (a/au), where a is the semimajor axis. 
7 is the mass ratio distribution parameter dN/dq oc e is the 
eccentricity distribution parameter dN/de oc (1 — e)^. 7 and e are 
calculated using maximum likelihood estimation. 


words, more compressive turbulence promotes more circu¬ 
lar orbits and more closely matched companion masses - 
but this is a very weak result. However, we can examine the 
combined distributions of semimajor axis a, mass ratio q 
and eccentricity e. Fig. shows the distribution of a, q and 
e alongside analytic fits. The parameters of these fits are 
given in Table These distributions include all the orbits 
of hierarchical systems. 

The semimajor axes range from 0.2 au to 2000 au. The 
lower limit is set by the resolution of the simulations (i.e. the 
radius of a sink particle). The upper limit is compatible with 
observations of young embedded populations (e.g. Taurus, 
Ophiuchus, etc, [King et al.||20I2b|a ). 

If we fit the distribution of mass ratios with a power 
law of the form dN/dq oc q^, we obtain 7 = 1.6 ± 0.2, 
implying a strong preference for companions of compara¬ 
ble mass. Young embedded populations also show a prefer¬ 
ence for companions of comparable mass, but it is somewhat 
weaker, viz. 0.2 ^7^1 ( Duchene V Kraus|20I3 ). To estab¬ 
lish whether the observations could be fit more closely with 
a specific value of (5 sol would require a much larger ensemble 
of simulations. 

If we fit the distribution of eccentricities with a reverse 
power law of the form dN/de oc (I — e)^, we obtain e = 
1.6 ± 0 . 2 , implying a preference for low-eccentricity, similar 
to the preferentially circular orbits of held M-dwarfs (see 
Duchene V Kraus||20I3 Fig. 4). 


5 DISCUSSION AND CONCLUSIONS 
5.1 Limitations 

First, we stress that the work reported here is a numeri¬ 
cal experiment, not an attempt to capture all the processes 
that ooccur in real prestellar cores. In particular we have 
simulated cores with a single mass, single radius, single non- 
thermal velocity dispersion, and single density prohle. The 
only parameters varied are the fraction of turbulent energy 
that is solenoidal, (5 sol, and the random seed that delivers 
different realisations, Xseed- Therefore the results cannot be 
applied to ah prestellar cores. Indeed, the results reported 
by LWHI4 suggest that cores with masses <C I M© tend to 
form single stars, whilst cores with > I M© tend to form 
discs and hlaments which fragment into multiple stars. The 
results presented here relate to the latter case. 

Second, we have not included magnetic helds in the sim¬ 
ulations. The three dimensional structure of a core’s mag¬ 
netic held is extremely difficult to extract from observations, 
and the evolution of the magnetic field requires the treat- 
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Figure 8. Histograms showing the distribution of orbital prop¬ 
erties integrated over all simulations (i.e. all (5 sol and Xseed)- The 
top panel shows the distribution of log^Q(a/au), where a is the 
semimajor axis. The middle panel shows the distribution of mass 
ratios, q. The bottom panel shows the distribution of eccentrici¬ 
ties, e. The red dashed lines are analytical fits to the data; fitting 
parameters are given in Table 


merit of non-ideal MHD effects, along with the detailed ion¬ 
isation chemistry of the matter. It is therefore beyond the 
scope of this investigation. 


tion; rather, turbulent energy is injected on many scales, in¬ 
cluding some that are smaller than, or comparable with, the 
core scale. In particular, when (if) a core condenses directly 
out of a shell swept up by an expanding nebula (Hii region, 
stellar-wind bubble or supernova remnant), there may be a 
predominance of compressive turbulent energy, a very small 
inertial range, and hence a low value of (5 sol, not so much 
disc fragmentation, and not so many low-mass stars. 


5.3 Comparison with similar work 


Girichidis et al. (2011) present adaptive mesh refinement 


simulations on scales intermediate between the molecular 
clouds simulated by, for example, [Bonnell, Clark Bate] 
(2008), and the cores treated in this work. They find that 


star formation occurs 25% earlier (with respect to the start 
time of the simulation) when the velocity field is purely com¬ 
pressive, as opposed to solenoidal. They do not report an 
increase in the number of low mass stars or brown dwarfs 
formed when the turbulence is purely solenoidal. This may 
be because their minimum resolvable length scale is tens 
of an, which is insufficient to capture the dynamics of disc 
fragmentation. 


5.4 Summary 

We have shown that the collapse and fragmentation of a 
core is influenced by the fraction of turbulent energy that 
is solenoidal, (5 sol, and hence - by default - by the fraction 
that is compressive. Specifically, as (5 sol is increased from 0 
(purely compressive) to 1 (purely solenoidal), the proportion 
of stars that form by filament fragmentation decreases and 
the proportion that form by disc fragmentation increases; 
at the same time the number of stars formed increases, and 
their mean mass decreases. The formation of massive cir- 
cumstellar discs requires (5 sol >1/3. With the limited num¬ 
ber of simulations that we have performed, we have been 
unable to establish any robust systematic trends in the mul¬ 
tiplicity statistics that derive from varying (5 sol. 


5.2 Implications for star formation 


Ssoh can only be determined observationally if one has the 
detailed three-dimensional velocity field, or if one has the 
detailed one-dimensional velocity field and a model, for ex¬ 
ample, assumed statistical isotropy of the velocity field (see 
Brunt Federrath|2014 ). On the scale of cores one does not 
have this information (because of instrumental limitations, 
selection effects, and confusion), and one cannot make these 
assumptions. 

The turbulent velocity field in a core is largely deter¬ 
mined by the flows that create it. If a core is created by tur¬ 
bulent fragmentation (e.g. Padoan Nordlund||20Q^ Hen- 


nebelle fc Chabrier[|2008[ |2009| and there is a large inertial 

range between the scale at which the turbulent energy is in¬ 
jected (say by galactic shear) and the scale of the core, Ssoh is 
likely to approach its thermal value, i.e. 2/3, and in this case 
disc fragmentation should be important. However, in reality 
there is probably not a clear division between the scale of 
energy injection and the inertial range leading to core forma¬ 
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